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We study the evolution of a population in a two-locus genotype space, in which the negative effects 
of two single mutations are overcompensated in a high fitness double mutant. We discuss how the 
interplay of finite population size, A'^, and sexual recombination at rate r affects the escape times 
tcsc to the double mutant. For small populations demographic noise generates massive fluctuations 
in tasc- The mean escape time varies non-monotonically with r, and grows exponentially as Infcsc ~ 
N{r — r*)^'''^ beyond a critical value r*. 



Point mutations at different loci of the genome af- 
fect the fitness of living organisms via complex intra- 
genomic correlations. These interactions, known as epis- 
tasis in population genetics, range from elementary pair- 
interactions to extended patterns reflecting the struc- 
ture of genetic networks Ij. An important form of 
epistasis termed reciprocal sign epistasis occurs when 
the deleterious effects of a single point mutation get 
(over)compensated by the beneficial effects of a sec- 
ondary mutation [21 [3] . The abundance of such elemen- 
tary motifs is expected from the nature of complemen- 
tary base pair binding in RNA |U [S] , but it also stands 
to reason that they act as building blocks of evolutionary 
processes in more complex 'fitness landscapes' [31 [5]. 

The effects of epistatic pair correlation can be conve- 
niently studied in a two-locus two-allele prototype sys- 
tem, i.e. a projection of the full genomic structure onto 
just two loci with alleles a, A and b, B, respectively. Here, 
the term 'locus' refers to a specific genomic position, and 
the 'alleles' a, . . . ,B stand for the nucleotides present at 
the loci. Throughout this paper we assume a fitness as- 
signment such that the double mutant AB is superior to 
the wildtype ab and both have higher fitness than the sin- 
gle mutants Ab, aB, whose fitnesses are moreover taken 
to be equal. When evolving under the joint influence 
of mutation, selection, sexual recombination and demo- 
graphic noise, this system exhibits dynamical phenomena 
on a variety of time scales. Of particular interest here 
is the 'switching' process occurring at very large time 
scales, when the effects of a deleterious single mutation 
at one of the two loci ab — > Ab or ab — > aB need to be 
overcome to reach a high-fitness doubly mutated config- 
uration {Ab/aB) AB. 

Several decades of research on the two-locus system 
notwithstanding [H |5l [7Hl6]. important aspects of the 
above compensatory mutation mechanism remain insuffi- 
ciently understood. While the switching process, termed 
stochastic tunneling in the literature [S] , is relatively well 
understood in asexuals [TTl [T2] , the behavior induced by 
sexual recombination is surprisingly complex. This is be- 
cause of the dual role played by recombination in this sys- 
tem. On the one hand, recombination of the unfavorable 



single mutants Ab + aB — >■ AB provides a channel of AB- 
generation that does not rely on mutation. On the other 
hand, this mechanism competes with back-recombination 
AB + ab ^ Ab, aB which breaks up ^i3-genotypes once 
they have formed. 

Detailed studies of the deterministic, infinite popula- 
tion dynamics show that the latter mechanism categori- 
cally wins in large populations. Beyond a critical value 
r* of the recombination rate r a stable stationary solu- 
tion localized at the wildtype genotype ab emerges |14j . 
and the escape to the double mutant AB is completely 
suppressed [13]. However, computer simulations of finite 
populations also display a parameter regime for r < r* 
where the escape is aided by recombination, such that an 
initial decrease of the escape time iosc eventually gives 
way to an increase at larger values of r [TO]. The effect 
of recombination on the speed of adaptation is thus seen 
to depend on the population parameters in a complex 
way, as has also been observed in studies using empirical 
fitness data [T7] . 

In addition to the recombination rate, the dynamics of 
the two-locus system is governed by the population size 
N, the mutation rate fi, and the strength of selection 
given by the typical scale of fitness differences s. Here 
we focus on populations that are moderately large, in 
the sense that NiJ,,Ns ^ 1, and subject to strong se- 
lection with s//i ^ 1 [TSj. Using a linearization of the 
full master equation near the initial wild type popula- 
tion, we derive an approximate expression for the full 
distribution of escape times for r < r* . This allows us 
to identify two fundamentally different switching scenar- 
ios. In populations that are smaller than a characteristic 
size to be specified below, the temporal bottleneck of the 
evolution is the appearance of the first few individuals 
of the AB populations, which is a rare event with an 
exponentially distributed waiting time. The most likely 
escape time is then much smaller than the typical time, 
and fluctuations in iosc are of the same order as the mean. 
In contrast, in large populations the evolution is limited 
by the growth of the AB population and the distribution 
of tcsc is sharply peaked around the mean. The non- 
monotonicity of ipse observed in [TUl is a feature of the 
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fluctuation-dominated regime which disappears in large 
populations. 

For r > r* the emergence of bistability in the deter- 
ministic dynamics implies that the problem becomes sim- 
ilar to the noise-driven escape from a metastable state. 
For this kind of problem path-integral methods akin to 
semi-classical quantum mechanics have recently been de- 
veloped [19], and we show that they can be applied in 
the present context as well. 

Model. We consider the dynamics of the two-locus 
system, as governed by the interplay of selection, mu- 
tation and recombination. The two low fitness genotypes 
Ah, aB are lumped into a single subpopulation. This is 
justified when the creation rate of single mutants sat- 
isfies N fjL ^ 1, such that the number of single mutant 
individuals is large compared to unity at all times and 
the recombination process Ah -\- aB — > (ab/AB) is not 
limited by the simultaneous presence of both parental 
types. The fitness of the three types is given by, re- 



spectively, ah -O- 1, (Ab/aB) o 1 — s^, AB o I + Sb, 
where Sd,b denote the selection coefficients of the dele- 
terious single mutants and the beneficial double mutant 
relative to the wild type ah, respectively. Mutation al- 
ters allelic content at a rate ii <C Sd, Sb- For instance, 

ah Ab/aB, where the factor of 2 accounts for the 
fact that the change of either allele, a or 6 generates a 
single mutant. For simplicity, we will neglect the effect 
of back mutation Ab/aB — >■ ah in our analysis of the 
escape time. In view of the smallness of the single mu- 
tant population (see below), this assumption is largely 
inconsequential. Finally, the random mating of individ- 
uals generates offspring by recombination at a rate r; for 
example, AB + ah (Ab/aB). 

The (Moran piT) master equation for the evolution of 
the three population sizes, no O ah, ni o Ah + aB, and 
ri2 O AB at constant total population = no + ni + 712 
reads as dtP{x,t) = H{x)P{x,t), where x — {xo,xi,X2), 
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[{E+ - l)d,(x) + [Er - l)6,,(x) + {E+Er^, - l)m,(x)] + V {E+EJ - l)r,,(x). 



(1) 



and the notation emphasizes the analogies to an imag- 
inary time 'Schrodinger equation'. In Q, = rii/N 
are the genotype 'frequencies', the 'operators' Ef f[xi) = 
e^v^'^i f(xi) — f{xi±N^^) act by translation by one in- 
dividual, and the rates di, b the death, birth, 
mutation, and recombination rates |21| . Equation ([T]) 
describes the evolution of the system at one individual 
death/birth event per time step. Alternatively, one may 
employ a Wright-Fisher approach where an entire gener- 
ation update is performed at each step . While this is 
numerically faster by a factor of 0{N), and has been used 
in parts of our computer simulations, the long rangedness 
of the Wright-Fisher equation in the variables (no, ni, n2) 
makes it difficult to handle analytically. 

Escape time. Figure [T] shows results for the time de- 
pendence of (no7"i7"2) obtained by simulation of the 
master equation ([T]). These profiles depend sensitively 
on the value of the parameter where the meaning of 

the rate ^ — 2^ — f-r|-^| will be discussed momentar- 
ily. For 'jN ^ 1 (upper panel) , the limiting factor for the 
escape process is the stochastic generation of a sufficient 
number of Ai?-individuals (cf. the inset). Once the AB- 
clone reaches a critical size, the fitness advantage triggers 
a fast sweep through the entire population. In the oppo- 
site case ''jN ^ 1, the quasi-deterministic increase of the 
a;2-population determines the escape time. 

We aim to describe the dynamical processes charac- 




FIG. 1. (Color online) Single runs of the escape process. Top 
(-yTV — 0.2): creation of the first AB-mutant determines es- 
cape time, bottom (7A'^ = 20): sweep through bulk of the 
population determines escape time. In the insets the pop- 
ulation size axes have been enlarged to highlight the initial 
stages of the process. 



terizing the respective temporal bottlenecks from the 
master equation. To this end, we first notice that 
starting from an initially un-mutated 'wild population' 
{xo,Xi,X2) = (1,0,0), a competition of mutation out of 
the a6-population and a counter-acting selection pressure 
Sd generates a population xi = xi = 2fi/sd at short time 
scales ^ Sd~^. In the parameter regime of interest here 
there are many individuals in a singly mutated state at 
any instance of time, Xi ^ N~^. During the evolution- 
ary stages determining the switching time, the frequen- 
cies cco — 1 and xi ~ xi do not change significantly 
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and the assumption of independence of P on (a;o,2;i) 
does not lead to qualitative errors. We thus substitute 
(xqjXi) ~ {l,xi) into ([T]) and set X2 = x = n/N for no- 
tational simplicity to obtain the linear master equation 
dtP{n,t) = H{n)P{n,t), where 

H{n) = (E- - l){R+n + 7) + {£+ - l)R^n (2) 

and R+ = ^(1 + Sfc + r/2), i?_ = ^(1 + 3r/2). The 
coefficient 7 has been defined above, and can now be 
understood as the effective rate of generation of double 
mutants, including contributions from both the muta- 
tional (~ ^ixi) and recombinational channel (~ rx\). 
The usage of the discrete variable n e N in ([2]) accounts 
for the importance of the 'quantization' of individuals to 
our further analysis. Starting from zero >li3-individuals, 
P{n,t = 0) = (5,1,0, the time scale for the generation of 
an A_B-clone is set by 7""'^. 

To obtain P{n,t) we adapt the general solution [22] 
for the generating function G{z,t) = X^^o -^"-^('^' ^) 
linear master equations to Eq. This yields G{z,t) = 
(1 + K(t)(l - z))-^, where A = "//R+ ~ N-/, K(t) = 
fli(gA_Rt _ g^^^ AR = R+-R_ > 0. Computing the 
inverse transform, we obtain 



P{n,t) = 



^ ' 3=0 
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For fixed n, P{n, t — >■ 00) scales to zero showing that the 
distribution flows towards large values of n. Writing the 
master equation in the form dtP = Jn-i — Jm we identify 
the current J„(t) = {R+n + "i)P{n,t) - R-{n + l)P{n + 
l,t) ~ ARnP{n,t) at which probability flows through a 
fixed reference value n 3> 1. The flow of the distribu- 
tion implies J dt Jn(t) = 1, which shows that J„(i) may 
be interpreted as the distribution /„(i) of escape times 
through n. We here consider n ~ N/2 as a reference value 
where the AB-population is about to take over, while the 
linearized approximation of the master equation has not 
yet become fully invalid. (For most parameter values, 
the ensuing threshold times are comparable to the times 
of full fixation, n = N.) For large n the product in ([3| 
can be expressed in terms of the F-function, which yields 
the result 



fnit) « AR 



K{tr 
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^F-i(A). 
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Two time scales can be extracted from Q. For t — > 00 
the distribution decays exponentially on the time scale 
ttaii = The second time scale is the time at which 

Q is maximal, given by iniax = (AR)^^ ln(Ai?n/7 -f 1). 
This is also the time when the solution of the rate equa- 
tion h = ARn + 7 corresponding to ([2| reaches the 
value n. Apart from the logarithmic factor we see that 
^max/itaii = A, showing that the parameter A distin- 
guishes the two dynamic regimes described above. For 



A <C 1 the distribution Q becomes purely exponential, 
and fluctuations in tosc are of the same order as the mean 
ttaii- In contrast, for A ^ 1 the distribution is sharply 
peaked around the most likely value tmaxj with fluctua- 
tions of order At with At/imax ~ N~^/'^. 

Inserting the explicit expressions for AR and 7 
it is straightforward to show that itaii displays a min- 
imum as a function of r provided Sd < Sf,/2, i.e. in 
cases where the fitness valley separating the genotypes 
ab and AB is relatively shallow. The minimum is lo- 
cated at rmin = Sb/2 — S(i < Sb/2. The maximal speedup 
due to recombination can be substantial, and is given by 
itaii(?-mi„)Ataii(r = 0) « ^ for ^ < s<j < Sfc 23j. By 
contrast, tmax increases monotonically with r. 

Metastability. Both time scales diverge when AR ^ 
{sb — t) vanishes as r — > s^. Within the simplified model 
defined by ([2]) , this reflects the emergence of a stable sta- 
tionary distribution centered around the wildtype geno- 
type ab, which has been found in previous studies of the 
deterministic model [HITnillJ. Figure [2] shows the mean 
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(triangles) along with the escape times ob- 



tained from the simulation of the Moran master equation 
([!]) (bars) , and from the solution of the corresponding de- 
terministic rate equations (dots). All three curves show 
a rapid increase of the escape time at the critical value 
r* ~ Sb- The rationale behind this behavior is that at 
strong recombination the 'reshuffling' AB + ab Ab, bA, 
outperforms the fltness advantage of the AB population, 
thus preventing its growth. 

We sketch the analysis of the model near criticality, 
r ~ r*. To this end we consider a Fokker-Planck ap- 
proximation to ([T]), i.e. a second order expansion in the 
'momentum' variables pi — dx^ ■ This results in an op- 
erator of the generic form H = —j^piFi + ^^^piXijpj, 
where the Fi denote the components of the deterministic 





FIG. 2. (Color online) Mean escape times as a function of r 
for Si, = 0.1, Sd = 0.025, = lO-*. Upper panel: N = 10^ 
lower panel. A'' = 10^. Error bars in lower panel illustrate the 
large fluctuations in the escape time, and inset highlights the 
minimum in tesc at r = rmin. 
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FIG. 3. (Color online) Main figure shows the logarithm of the 
escape time, normalized to its value at r = r*, as a function of 
r — r* in double-logarithmic scales. The escape time was de- 
fined as the time when the AB-clone reaches 80% of the total 
population. Black dots correspond to simulations of a popu- 
lation of size A'^ — 2^^ ~ 6.7 x lO'^ with parameters Sb — 0.1, 
Sd = 0.025, and /i = 10~^. All data are averaged over 10"" 
runs. The straight line has slope 3/2. Inset illustrates the de- 
terministic one-dimensional dynamics of the ^_B-population. 
The vertical graph on the left shows the shape of the func- 
tion F{x) for r < r* (dashed line) and r > r* (dotted line). 
The horizontal graph shows the evolution of the fixed point 
structure with r. For r < r* the only fixed point is near 
X2 = X ^ 1; for r > r* another pair of fixed points (Xs,a;*) 
emerges near a; ~ 0. 



drift term and Xij is the diffusion matrix. Eliminat- 
ing one frequency by normalization, xo = 1 — a;i — X2^ 
the problem becomes two-dimensional with i = 1,2. 
The quickly equilibrating frequency Xi ~ — of the sin- 
gle mutant population follows from stationarity under 
drift, Fi{xq,xi) — 0, which gets us to the effective one- 
dimensional Hamiltonian, H = —j^pF+ i^j^Xp^^ where 
we have set X2 = x and P2 ^ P for notational simplicity. 
Then F{x) = F2{x,Xi), X{x) = X22{x,Xi) ~ x{l — x) 
describe the effective dynamics of the ^_B-population. 

The singularities in the escape times reflect the emer- 
gence of a pair of stable {x*) and unstable (a;*) fixed 
points, i^(a;*_„) = 0, < x* < a;* < 1, close to the 
boundary a: = (inset of Figjs]). The unstable fixed point 
X* defines a 'recombination barrier' [4] which blocks the 
fixation of the fitter population out of a small number of 
initial individuals x(0) < a;* . Based on general results for 
the escape from metastable states obtained within a large 
deviations ('WKB') approach [IS], one expects the escape 
time to grow exponentially with N as tesc ^ exp[CiV5^], 
where C > is a constant and (5 ~ x* — x* . The analysis 
of the deterministic Wright-Fisher version of the prob- 
lem shows that 5 ^ [r — r*)^/^ [T3], and we predict 
that Use - exp[CiV(r - r*)3/2]. Figure [s] shows data for 
the escape time obtained by simulating a Wright-Fisher 
process [17] at values of r slightly above r*. The value 



r* ~ 0.09172 was calculated using expressions derived in 
[l4] for the deterministic (iV — )■ oo) limit. The results 
are seen to be in good agreement with the prediction 
lnie.c^(r-r*)3/2. 

Conclusion. To summarize, we have presented an anal- 
ysis of a paradigmatic two-locus model of population ge- 
netics. We have seen how, depending on the population 
size, recombination may speed up or delay the evolu- 
tion towards the high fitness state. A challenge for fu- 
ture work is to develop the WKB-type approximation 
for r > r* into a fully quantitative theory, which can be 
used to predict the size of the recombination barrier in a 
specific biological setting. 
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